function [K11,K12,K22, b1,b2] = elastic_T3(V, T, Quad_Order, FE_Type, FE_Order, f, Young, nu,  varargin)
% function [K11,K12,K22, b1,b2] = elastic_T3(V, T, Quad_Order, FE_Type, FE_Order, f, Young, nu, varargin)
%  ?��Ԫ��ɢ
% T3 assume that the geom mapping is linear
%
% The Mesh is Given in (V,T)
% Quad_Order is the precision of the quadrature
% FE_Type and FEO_Order is the informaiton for FE
% lambda, mu, f, g are informations of the planar elastic problems
%
% The returned value are ElementMatrix and ELement Right Hand Side 
% which are stored elementwisely and block-structured!
%

lambda = Young*nu/((1+nu)*(1-2*nu)); 
mu = .5*Young/(1+nu);

 n_elem = size(T, 1);
 n_basis = (FE_Order + 1)*(FE_Order + 2)/2;
 
% reserve spaces for the return values
%  ElemMatrix = zeros(n_basis, n_bais, n_elem);
%  ElemRHS = zeros(n_bais, n_elem);
 K11 = zeros(n_basis, n_basis, n_elem); 
 K12 = zeros(n_basis, n_basis, n_elem);
 K22 = zeros(n_basis, n_basis, n_elem);
 b1 = zeros(n_basis,n_elem);  
 b2 = zeros(n_basis,n_elem);
 
 %%%%%
 % 1. the quadrature informations
 [qw, qp] = quad_rule(Quad_Order);
 
 %%%%%
 % 2.  the shape function value and their first order gradients
 [phi, dphidxi, dphideta, ddxi, ddxieta, ddeta] = fe_basis(FE_Type, FE_Order, qp(:,1),qp(:,2));
 
 %%%%%
 % 3. the physical coordinates for all quadrature points
 r = qp(:,1);   s = qp(:,2);   t = 1 - r - s;
 px = (r*V(T(:,1),1)' + s*V(T(:,2),1)' + t*V(T(:,3),1)');
 py = (r*V(T(:,1),2)' + s*V(T(:,2),2)' + t*V(T(:,3),2)');

 [f1,f2] = feval(f, px, py, Young, nu, varargin{:});
  
 %%%%%
 % calculate the Jacobian and derivatives on each elements
 % dphidx, dphidy sized n_q_pts - by - n_basis x n_elem
%  n_q_pts = length(qw);
 xxi   = V(T(:,1),1) - V(T(:,3),1);
 xeta = V(T(:,2),1) - V(T(:,3),1);
 yxi   = V(T(:,1),2) - V(T(:,3),2);
 yeta = V(T(:,2),2) - V(T(:,3),2);
 jacobian = xxi.*yeta - xeta.*yxi;
 a11 = yeta./jacobian;  a12 = - yxi./jacobian;
 a21 = -xeta./jacobian; a22 = xxi./jacobian;
 
 %%%%%
 % 4.  assebmle the planar elastic stiff matrix element-wise
 for k = 1:n_elem
     weight = 0.5*jacobian(k)*qw*ones(1,n_basis);
     
     phiwei = weight.*phi;
     dphidx = a11(k)*dphidxi + a12(k)*dphideta;
     phixwei = weight.*dphidx;
     dphidy = a21(k)*dphidxi + a22(k)*dphideta;
     phiywei = weight.*dphidy;

     K11(:,:,k) = (lambda + 2*mu)*dphidx'*phixwei + mu*dphidy'*phiywei;
     K12(:,:,k) = (lambda +    mu)*dphidx'*phiywei;
     K22(:,:,k) = (lambda + 2*mu)*dphidy'*phiywei + mu*dphidx'*phixwei;
     
     b1(:,k) = phiwei'*f1(:,k);
     b2(:,k) = phiwei'*f2(:,k);
 end
 
end